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ABSTRACT 


A  fully  Sinc-Galerkin  method  for  recovering  the  spatially  varying  stiffness  and  damp¬ 
ing  parameters  in  Euler-Bernoulli  beam  models  is  presented.  The  forward  problems  are 
discretized  with  a  sine  basis  in  both  the  spatial  and  temporal  domains  thus  yielding  an  ap¬ 
proximate  solution  which  converges  exponentially  and  is  valid  on  the  infinite  time  interval. 
Hence  the  method  avoids  the  time-stepping  which  is  characteristic  of  many  of  the  forward 
schemes  which  are  employed  in  parameter  recovery  algorithms.  Tikhonov  regularization  is 
used  to  stabilize  the  resulting  inverse  problem,  and  the  L- curve  method  for  determining 
an  appropriate  value  of  the  regularization  parameter  is  briefly  discussed.  Numerical  exam¬ 
ples  are  given  which  demonstrate  the  applicability  of  the  method  for  both  individual  and 
simultaneous  recovery  of  the  material  parameters. 

1This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NAS1-18605  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  2366S. 
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1  Introduction 


In  the  modeling  and  control  of  large  flexible  structures,  one  is  often  required  to  numerically 
recover  one  or  more  material  parameters  given  data  measurements  at  various  points.  Al¬ 
though  these  structures  are  in  general  very  complex,  in  many  cases  the  essential  features  can 
be  developed  by  considering  a  fixed  Euler-Bernoulli  beam  which  is  assumed  to  have  Kelvin- 
Voigt  damping.  In  this  paper,  a  fully  Sinc-Galerkin  method  is  presented  for  the  numerical 
recovery  of  the  stiffness  parameter  El  and  the  damping  parameter  cpl  in  the  state  space 
model 


C(EI,cpI)u  =  f(x,t),  0  <  x  <  1,  t  >  0 


ti(0,  <)  =  tt(l,  t)  =  0,  f  >  0 

§>0 -£(1.0-0.  <>0 

r\ 

u(x,0)  =  0)  =  0,  0  <  x  <  1 


(1.1) 


with 


given  measurements  of  the  data  at  the  points  {(®p,  in  the  domain  (0,1)  x  2R+. 

From  physical  considerations,  it  is  reasonable  to  let  the  admissible  parameter  set  Q  be  defined 

by  a 

Q  =  j(£/,cD/)  €  ]I  #  2(°.  !) :  EI(X)  >EIo>0,  cDI(x)  >  o| 

(see  [5]).  With  this  definition,  the  existence  of  a  unique  solution  u  to  the  forward  problem 
can  be  obtained  on  any  fixed  time  interval  [0,  r],  r  >  0,  for  /  sufficiently  smooth. 

The  “idealized”  parameter  recovery  problem  can  then  be  formulated  as  follows:  determine 
q  =  (El,  cpl)  €  Q  such  that 

Cu(;;q)  =  d  (1.2) 

where  u(-,  • ,q )  =  C~1(q)f  denotes  the  parameter-dependent  state  solution  and  d  is  the  data. 
The  observation  operator  C  is  given  by 

(i.3) 


i 


and  hence  Cu  represents  point  evaluations  of  the  solution.  Note  that  (1.2)  can  be  written  as 
the  operator  equation 

K(q)  =  d  (1.4) 

with  the  nonlinear  operator  K  given  by 

K(q)  =  CC~\q)f  . 

For  reasons  similar  to  those  discussed  in  [11],  the  problem  (1.4)  is  ill-posed  in  the  sense 
that  solutions  q  (provided  they  even  exist)  may  not  depend  continuously  on  the  data  d. 
Consequently,  some  sort  of  regularization  (i.e.,  stabilization)  is  required  to  obtain  an  accurate 
approximation  for  q. 

The  regularization  technique  that  is  used  is  Tikhonov  regularization  [13],  and  the  problem 
(1.4)  is  replaced  by  the  minimization  problem 

min  Ta(q)  (1.5) 

where 

Here  a  >  0  is  a  regularization  parameter  which  controls  the  tradeoff  between  goodness  of 
fit  to  the  data  and  stability.  The  penalty  functional  J ( q )  provides  stability  and  allows  the 
inclusion  of  a  priori  information  about  the  true  parameters. 

Due  to  the  infinite  dimensionality  both  of  Q  and  of  the  state  space,  the  problem  (1.5)  is  an 
infinite  dimensional  minimization  problem.  In  order  to  develop  a  practical  numerical  scheme, 
the  problem  must  be  replaced  by  a  sequence  of  finite  dimensional  problems;  that  is,  one 
must  approximate  the  operator  K  and  minimize  the  functional  Ta  over  a  finite  dimensional 
admissible  subspace  of  Q. 

The  evaluation  of  K(q)  requires  the  solution  of  (1.1).  Similar  partial  differential  equa¬ 
tions  must  be  solved  to  obtain  the  components  of  the  derivative  fC'(q).  The  construction 
of  an  approximate  solution  to  these  forward  problems  commonly  begins  with  a  Galerkin 
discretization  of  the  spatial  variable  with  time-dependent  coefficients.  This  yields  a  system 
of  ordinary  differential  equations  which  is  solved  via  differencing  techniques.  For  problems 
with  nontrivial  cp/,  it  is  noted  in  [1]  that  the  equations  are  moderately  stiff  and  routines  for 
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stiff  systems  must  be  employed,  thus  adding  to  the  expense  of  the  algorithms.  This  difficulty 
is  further  augmented  by  the  fact  that  the  time-stepping  must  be  repeated  at  each  step  in  the 
minimization  of  (1.5).  A  final  difficulty  lies  in  the  need  to  interpolate  at  data  points  which 
do  not  coincide  with  the  nodes  of  the  ODE  solver. 

In  contrast,  the  method  of  this  work  implements  a  Galerkin  scheme  in  time  as  well  as 
space,  thus  bypassing  many  of  the  difficulties  associated  with  time-stepping  methods  in  the 
context  of  inverse  problems.  Corresponding  results  for  the  heat  equation  can  be  found  in 
[7],  and  a  detailed  discussion  of  fourth-order  results  involving  the  recovery  of  El  in  models 
with  no  damping  is  given  in  [11]. 

The  fully  Sinc-Galerkin  method  in  space  and  time  has  many  salient  features  due  both 
to  the  properties  of  the  basis  functions  and  the  manner  in  which  the  problem  is  discretized. 
Perhaps  the  most  distinctive  feature  of  the  method  is  the  exponential  convergence  rate 
when  solving  the  corresponding  forward  problems.  Furthermore,  the  judicious  choice  of  a 
conformal  map  provides  approximate  solutions  which  are  valid  on  the  infinite  time  interval 
rather  than  only  on  a  truncated  time  domain.  Finally,  the  discrete  system  requires  no 
numerical  integrations  to  fill  either  the  coefficient  matrices  or  the  right-hand  side  matrix. 
All  three  features  prove  to  be  advantageous  when  solving  the  forward  problems  and  hence 
the  inverse  problem. 

In  Section  2,  the  Sinc-Galerkin  system  for  the  forward  problem  is  considered  and  imple¬ 
mentation  details  are  discussed.  The  forward  results  are  then  incorporated  into  the  finite 
dimensional  minimization  problems  in  Section  3  with  the  discussion  centering  around  the  con¬ 
struction  of  the  Tikhonov  functional.  Numerical  results  are  presented  in  Section  4  along  with 
a  brief  outline  of  the  “L-curve”  technique  for  determining  the  regularization  parameter  a. 
Examples  are  given  which  demonstrate  the  recovery  of  the  individual  parameters  EJ  and 
cp/  as  well  as  the  simultaneous  recovery  of  both  parameters. 
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2  The  Forward  Problem 

Consider  the  forward  problem 


Cu(x,  t)  =  /(x,  f),  0<x<l,  t>0 

u(0,  f)  =  u(l,t)  =  0,  i>0 

|(o, 0  =  ^(1, 0  =  0,  i>0 

u(x}0)  =  -~(x,0)  =  0,  0  <  x  <  1 

„  /  x  diu  /„r,  .3Ju.  T,  .  d3u  . 

Cu(x,  t)  =  — (x,  0  +  gji  o  + 


Since  a  thorough  derivation  of  the  Sinc-Galerkin  method  for  problems  of  this  type  (with 
cdI  =  0)  is  given  in  [10],  the  following  discussion  contains  only  that  material  which  is 
needed  for  the  construction  of  the  associated  matrix  system. 

A  fully  Sinc-Galerkin  method  for  the  approximation  of  the  solution  of  (2.1)  may  be  briefly 
summarized  as  follows.  For  <f>(x)  =  In  T(u;)  =  ln(w)  and  positive  hx  and  ht  ,  define 

Si(x)  =  S(i,  hx)  o  ^(x)  =  jinc  _  (2-2) 

$(<)  =  S(},  h,)  o  T(t)  =  sine  (T(‘\yh')  (2  3) 


where 


.  .  sin(xx) 

atnc(x)  s  — ' — — oo  <  ®  <  oo  . 

7TX 


The  basis  is  then  taken  to  be  {5^}  with 

Sij(x,t)  =  5j(x)Sj(t), 

and  the  approximate  solution  is  defined  by  way  of  the  tensor  product  expansion 


N,  Nt  m,  =  Mx  +  N,  +  1 

Wm.mi(3!i0  =  S  (2-5) 

mt  =  Mt  +  Nt- fl. 
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The  m,  •  mt  unknown  coefficients  {«»,■}  are  determined  by  orthogonalizing  the  residual  with 
respect  to  the  set  of  sine  functions  { SpS J M,  •  This  yields  the  discrete  Galerkin 
system 

(£u„.m,-/,sps;)  =  o  (2.6) 

for  p  =  —Mx,  •••  ,NX  and  q  =  —  Mt,  •  •  • , Nt.  The  inner  product  (•,  •)  is  taken  to  be 

(F,G)  =  J  J  F(xyi)G(x,t)w(xtt)dxdt  (2-7) 

where  the  weight  is 

w(xyt)  =  w(x)w*(t)  =  (^'(x))'*(T (t))~*  .  (2.8) 

A  thorough  discussion  motivating  this  choice  of  weight  can  be  found  in  [10]. 

The  expressions  (2.1),  (2.7),  and  (2.8)  are  then  combined  to  form  the  system 

f°°  f 1  d 3 

Jo  Jo  Q^(Um>m‘)^pSlww* dxdt 

+C  /„'  &  ■)) 

(2.9) 

f°°  /»  d3  /  d3  \ 

+/„  Jo  a^{CD'wm{u — ))  s^urw-ixdt 

—  J  J  fSpSgWw'dxdt 

for  p  =  -Mx,  •  •  • ,  Nx  and  q  =  —Mty  •  •  • ,  Nt. 

In  anticipation  of  the  parameter  identification  problems  which  motivate  this  analysis,  the 
terms  El  and  cqI  in  (2.9)  are  expanded  as  linear  combinations  of  weighted  sine  functions 
with  four  Hermite-like  algebraic  terms  added  to  accommodate  the  potentially  nonzero  func¬ 
tion  and  derivative  values  of  El  and  cqI  at  x  =  0  and  x  —  1.  Specifically,  this  parameter 

basis  is  taken  to  be  {ipk}k=-Mm  w^h 

* 

x(l  —  x)a,  k  =  —Mx 

(1  -  ®)a(2x  +  1),  fc=  -Mx  4- 1 

M*)  =  '  vE(x)Sk(x),  k  =  -Mx  -f  2,  •  •  • ,  Nx  -  2  (2.10) 

®a(2(l  -  *)  4- 1),  k  =  Nx-  1 
-(1  -  x)xa,  k  =  Nx  . 
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Here  S*(x)  =  S(k,hx)  o  <f>(x)  and  the  basis  weight  is  taken  to  be 


vb(x)  =  w(x)  =  [®(1  —  x)]i  . 


(2.11) 


The  finite  dimensional  approximations  for  El  and  cDI  are  then  given  by  the  expansions 


and 


EIm.{x)  =  53 

fc=-M. 

N. 


(2.12) 


(2.13) 


CD^m.(x)  =  53  Ckrl>k{x)  . 

k=-Mm 

In  the  forward  problem,  the  coefficients  {e*}  and  {c*}  are  known  whereas  in  the  corre¬ 
sponding  parameter  recovery  problems,  they  are  unknown  and  are  determined  via  methods 
to  be  discussed  in  Section  3.  The  number  of  basis  functions  used  in  the  expansions  is  chosen 
so  as  to  guarantee  a  square  spatial  coefficient  matrix. 

A  quick  note  should  be  made  concerning  the  choice  of  parameter  basis  and  the  manner  of 
expanding  EImm  and  co/m..  The  two  derivative-interpolating  boundary  basis  functions  are 
added  so  that  these  expansions  are  the  same  as  those  used  with  cantilever  or  free  boundary 
conditions.  The  choice  of  (2.11)  for  the  basis  weight  is  certainly  sufficient  and  proves  to  be 
beneficial  when  incorporating  this  forward  scheme  into  a  numerical  method  for  solving  the 
parameter  recovery  problem  as  described  in  Section  3. 

Sine  quadrature  is  used  to  evaluate  the  integrals  in  (2.9)  and  hence  derive  a  discrete 
system.  For  details  of  the  quadrature  rule  and  conditions  governing  its  error  bound,  see 
[12];  for  the  purposes  of  this  work,  however,  it  suffices  to  state  the  sine  quadrature  results  as 
follows.  Let  T  be  (0,1)  or  (0,  oo)  when  x  =  <t>  or  T,  respectively.  If  F  is  analytic  on  T  and 
suitably  bounded,  and  if  there  exist  positive  constants  K,a,  and  (3  such  that 


F(r) 


X'(T) 


<K( 


t  e  ^((-oo^o)) 

e"^lx(T)',  tG^>([0,  oo)) 


(2.14) 


where  V*  —  X-1»  then  for  h  >  0  sufficiently  small 


Lw-hkM 


j=~M  X'(zj) 

The  sine  gridpoints  are  given  by  Zj  =  ip(jh)  =  x_1(i^) 


<  +  -e-"k  + 

a  P 


(2.15) 
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The  expansions  (2.12)  and  (2.13)  are  substituted  into  (2.9),  and  integration  by  parts  is 
used  to  transfer  the  derivatives  onto  the  product  SpwS*w*.  As  detailed  in  [10],  the  weight 
choice  (2.8)  guarantees  that  all  boundary  terms  vanish.  The  resulting  integrals  are  then 
evaluated  via  (2.15).  For  Y  =  El  or  CpJ  and 

y(x)  =  Y(x)  —  V(0)[(1  —  x)J(2z  -f  1)]  —  y’(l)[®a(2(l  —  a)  +  1)] 


-5"(0)[x(l  -  *)»]  -  K'(l)[-(1  -  x)x»J  , 


the  requirement 

|.V(a:)it(z,  J)|  <  Kxa**(  1  —  x)fi+^{y+^e~Sl  (2.16) 

guarantees  the  decay  needed  to  truncate  the  infinite  quadrature  rule  as  specified  by  (2.14). 
With  a,  0, 7  and  8  specified  and  Mx  chosen,  the  choices 


and 

for  the  stepsizes  and  summation  limits  balance  the  asymptotic  errors  to  at  least  order 
0(e(~*daMm^).  This  rate  results  from  the  presence  of  a  sine  function  in  the  integral.  In 
the  above  expressions,  [•]  denotes  the  greatest  integer  function.  Note  that  the  +1  is  unnec¬ 
essary  when  jjMx  or  ~MX  are  integers. 

Given  MX)  Nx,  Mtj  Nt  and  h  =  hx  =  ht  as  defined  above,  the  discrete  system  for  (2.1)  is 


Ae:UCf  +  A,c,UAl  +  CJUAl  =  G. 


(2.17) 


and 


and 

[p] a  =  /(*.•>  tj)  ■ 

It  should  be  noted  that  the  ordering  of  the  coefficients  t in  U  mimics  that  used  in  most 
standard  time-differencing  schemes.  This  is  a  matter  of  convenience  since  the  Sinc-Galerkin 
method  is  not  bound  by  any  specific  ordering  of  the  grid. 

As  shown  in  [8],  the  mt  x  mt  temporal  matrices  An  and  At2  are  given  by 


=  [h,m  ~  ?'] 

where  I  denotes  the  mt  x  rrij  identity,  and  JW  and  are  defined  componentwise  by 

0,  j  =  q 


3  ’ 

(_2)(-l)j-9 

O'-?)3  ’ 


3  =  ? 
^  ?  • 


The  m,  x  m*  spatial  matrices  Agj  and  AeDj  have  the  form 


Ab,  =  +  2#<»X>(?„o)  +  . 

dI  =  [$<’>©(  i>„„)  +  2f<’)p(p,c.,)  + 
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where 


P*(«>  =  t  =  0,1,2  , 

PU*)  =  /=  0,1,2 

with  c  =  [c_M.,  •  •  •  ,c^.]r  and  c  =  [c_w.,  ■  •  • ,  c^.]T.  Finally,  the  matrices  j  =  2,3,4 
and  tfM,  l  =  0, 1, 2  are  defined  componentwise  by 

l*WU  =  -^(S,«0W(*i)  (2.18) 

and 

[¥W]i»  =  #(*.')  (2.19) 

with  the  notation  on  the  right-hand  sides  of  (2.18)  and  (2.19)  indicating  the  j-th  and  /- th 
derivatives,  respectively.  The  expansions  of  and  tyW  in  terms  of  more  fundamental 
matrices  can  be  found  in  [10]. 

For  given  El  and  cpl,  one  then  needs  to  solve  the  linear  system  (2.17)  for  the  matrix  U. 
As  shown  on  page  414  of  [6],  (2.17)  is  algebraically  equivalent  to 

Au  =  {Ct  ®  Aei  4-  Ati  ®  ACt)i  -f  At2  ®  Ct}  co{U)=co{G)  (2.20) 

where  the  vector  u  =co(U)  is  the  concatenation  of  the  mx  x  m,  matrix  U  obtained  by 
successively  “stacking”  the  columns  of  U,  one  upon  another,  to  obtain  an  (mx  •  mt)  x  1 
vector.  The  system  (2.20)  can  be  solved  directly  via  any  of  the  decomposition  methods  that 
are  available  for  large  linear  systems.  It  should  be  noted  however,  that  for  large  values  of  mx 
and  mg,  the  ill-conditioning  of  the  (mx  •  mt )  x  (mx  •  mt)  matrix  A  must  be  considered  when 
devising  schemes  for  solving  the  system  (2.20).  One  facet  of  current  research  is  directed  at 
devising  linear  algebra  techniques  which  will  facilitate  the  solution  of  this  system. 
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3  The  Finite  Dimensional  Minimization  Problem 


As  noted  in  the  introduction,  the  minimization  problem 

where 

Ui)  =  j{||JC(?)  -  df  + 

is  infinite  dimensional  and  thus  must  be  replaced  by  a  sequence  of  finite  dimensional  problems 
before  a  viable  numerical  scheme  can  be  developed.  Following  from  (2.12)  and  (2.13),  the 
approximating  admissible  parameter  sets  are  taken  to  be 

Qm,  =  {(5fm,/ofm,)  •  EImm(x)  —  Y^k=-l \fm  cklpk(x), 

CDIm.(x)  =  hMX)} 

with  the  basis  {ipk}  defined  in  (2.10).  For  qmm  =  (EImm,coImm),  the  associated  finite  dimen¬ 
sional  optimization  problem  can  then  be  formulated  as 

*3 1.  *•(*"•)  <31> 

where 

t.( ?m.)  =  j  {||  *(»„.)  -  rfll'  +  .  (3.2) 

Note  that  in  solving  the  minimization  problem  (3.1),  one  is  actually  solving  for  the  vec¬ 
tors  c  =  [c_ m.  ,  •  •  • ,  ctf.]  G  2Rm’ ,  c  =  (c_m.  ,  •  •  • ,  cat.]  G  Rmm ,  or  c  =  [c,  c]  G  depending 

on  whether  one  is  solving  for  EIma,  or  both  parameters  simultaneously. 

With  np  and  n,  specifying  the  number  of  spatial  and  temporal  observation  points,  respec- 

A 

tively,  the  approximation  K(qm „)  to  IC(q)  is  obtained  by  applying  the  observation  operator  C 
in  (1.3)  to  um.m,  in  (2.5).  If  the  set  of  observation  points  {(xp,  ^)}p=i’..  ’n*  can  be  represented 

A 

as  a  tensor  product  of  spatial  and  temporal  points,  then  K(qmm)  has  the  representation 

*(?«,.)  =  C  c^U)  (3.3) 

where  the  matrix  U  solves  (2.17).  C  is  an  (np  n,)  x  (m,  Tnt)  evaluation  matrix  which  can  be 
formulated  as  follows.  Define  the  np  x  m,  spatial  evaluation  matrix  Ex  to  have  components 

=  ^  V  ^  np>  —A/*  <  i  <  Nz 
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and  let  the  n,  x  mt  temporal  evaluation  matrix  Et  have  the  components 

[EtUj  =  S;(tq),  1  <q<nqt  -Mt  <  j  <  Nt. 

Then 

C  =  Et  ®  £*. 

It  is  noted  that  if  the  set  of  observation  points  is  not  rectangular  as  described  above,  then 
point  evaluation  can  be  done  directly  via  (2.5).  This  latter  option  is  less  efficient,  however, 
than  that  defined  in  (3.3). 

The  form  of  the  discrete  penalty  functional  J  depends  upon  whether  one  is  solving  for 
Elm.,  CDlmm,  or  both  simultaneously.  In  the  first  case,  the  discrete  penalty  functional  is 
taken  to  be 

A*n.)  =  JjECm(x)]2dx  +  e  JjEIm.(x)]2dx  «  ctQc 
where  the  m,  x  m*  matrix  Q  =  Qd  -f  Qj  has  components 

[Qa]kt  «  fl  Wxytftx)dst  ~Mt  <  k,e<Nt 

Jo 

and 

[Q/]m  «  /  il>k{x)ilu{x)dx,  -Mx  <k,£<  Nx. 

Jo 

The  exact  representations  for  the  matrices  Qj  and  Qf  are  given  in  [11].  Similarly,  if  EImm  is 
known  and  is  unknown,  an  appropriate  penalty  functional  is 

«J(9m.)  =  /  [co/m.(*)]J<ia:  +  t  I  [cD/m.(x)]adl  «  CT Qc 
Jo  Jo 

with  Q  defined  as  above.  Finally,  in  the  case  where  both  parameters  are  unknown,  the 
discrete  penalty  functional  can  be  taken  to  be 

*J(?m.)  »  ctQc 


where  the  2m*  x  2m*  matrix  Q  is  given  by 


Qf  +  Qd  0 
0  Qf  +  Qd 
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Although  the  matrix  Q  in  each  case  is  not  sparse,  as  shown  in  [11]  it  is  very  efficient  to 
construct  since  the  component  matrices  are  also  needed  in  the  forward  solver.  Moreover,  for 
e  >  0  the  matrix  Q  is  symmetric  and  positive  definite  and  hence  has  a  Cholesky  decompo¬ 
sition  Q  =  RTR  where  R  is  upper  triangular.  If  c  is  used  to  denote  c,  c,  or  c,  then  it  follows 
that  the  penalty  term  J(qma )  yields  the  quadratic  form 

ctRtRc=\\Rc\\2  (3.4) 

where  ||  •  ||  denotes  the  Euclidean  norm.  This  representation  for  the  penalty  functional  is 
particularly  useful  both  when  implementing  a  scheme  to  solve  the  minimization  problem  and 
when  plotting  the  “L-curve”  to  determine  a  suitable  regularization  parameter  a  (see  [3]). 

A 

To  highlight  the  dependence  of  the  functional  Ta  in  (3.2)  on  the  unknown  vector  c,  let 

K(c)  =  K(qmm)  =  C  co(U(c)) 

where  17(c)  solves  (2.17)  for  a  given  c.  Noting  (3.4),  the  optimization  problem  (3.1)  can  be 
replaced  by 

mm  Ta(c)  (3.5) 

where 

Uc)  =  |{||A-(5)  -  df  +  °||Rc||’}.  (3.6) 

To  obtain  a  minimizer  for  the  nonlinear  functional  Ta{ c),  a  quasi-Newton  trust  region  iter¬ 
ation  [2] 

Cfc+l  =  Cfc  +  Sfc 

is  used.  Here  solves  the  quadratic  programming  problem 

m.in  \  +  K'(e^3  ~  dW3  +  a\\Rfa  +  *)HJ} 

subject  to  |H|  <  6k  with  K\ck )  denoting  the  Jacobian  of  K  at  c*.  The  trust  region  radius  is 
chosen  so  that  Ta(c )  has  sufficient  decrease  at  each  iteration  to  guarantee  convergence  to  a 
local  minimizer  of  Ta  (for  further  details  about  the  theory  and  implementation  of  the  trust 
region  algorithm,  see  [2]  or  [4]). 
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An  important  numerical  issue  in  the  implementation  of  the  trust  region  scheme  is  the 
formulation  of  the  derivative  of  the  operator  K.  Here  the  derivative,  or  Jacobian,  is  the 
matrix  whose  i/-th  column  is  given  by 

m<0 ],  =  Km  ~[K(c  +  Te„)  -  K(c)\ 
where  the  standard  unit  vector  e„  has  components 


1, 

0, 


k  =  v 
k  ^  v  . 


In  the  examples  of  the  next  section,  the  Jacobians  are  calculated  with  a  standard  forward 
difference  scheme.  This  scheme  is  easy  to  implement  and  accurate  enough  for  the  purposes 
of  the  method.  If  further  efficiency  is  desired,  a  directional  derivative  scheme  such  as  that 
described  in  [7]  can  be  used. 


4  Implementation  and  Numerical  Examples 

The  three  examples  in  this  section  demonstrate  the  use  of  the  Sinc-Galerkin  method  for 
recovering  the  individual  parameters  El  and  cqI  as  well  as  the  simultaneous  recovery  of 
both  parameters.  In  each  case  the  state  solution  is  u(x,t)  =  ®(1  —  x)sin(4 xx)tae_‘  and  the 
true  material  parameters  are  EI(x)  =  1  +  x  +  x(l  —  x)  and  cpJ(x)  =  1  +  sin(xx).  For  these 
functions,  the  choices  a  =  fi  =  7  =  \  and  5=1  satisfy  the  decay  condition  (2.16).  In 
all  three  cases,  the  dynamics  of  the  problem  are  assumed  to  be  modeled  by  (1.1)  with  the 
forcing  function  /(®,  t )  consistent  with  the  state  solution  and  the  true  stiffness  and  damping 
parameters. 
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In  the  first  example,  El  is  considered  known  and  c^I  is  numerically  recovered  using  the 
methods  of  this  paper.  The  roles  are  reversed  in  the  second  example  and  El  is  numerically 
recovered  while  cjpJ  is  considered  known.  The  final  example  demonstrates  the  numerical 
recovery  of  both  unknown  parameters.  It  should  be  noted  that  a  comprehensive  set  of 
examples  demonstrating  the  application  of  the  Sinc-Galerkin  method  to  models  of  undamped 
beams  (cp/  ~  0)  can  be  found  in  [11]. 

A  very  important  practical  consideration  is  the  choice  of  the  regularization  parameter  a 
for  a  given  (error-contaminated)  data  set.  If  the  error  in  the  data  is  discrete  and  random, 
then  under  certain  conditions  the  method  of  Generalized  Cross  Validation  (GCV)  can  be 
used  to  determine  a  suitable  value  of  a  (see  [14]).  A  second  method  for  determining  the 
regularization  parameter  is  to  plot  the  norm  of  the  penalty  functional,  ||/?ca||,  versus  the 
norm  of  the  residual,  ||/f(ca)  —  d||  (see  [3]  or  [9]).  Here  ca  denotes  the  solution  to  (3.5).  In 
this  way,  one  can  qualitatively  get  an  idea  of  the  compromise  between  the  minimization  of 
these  two  quantities.  The  scheme  for  determining  the  “optimal”  regularization  parameter 
consists  of  finding  those  values  of  a  such  that  (||Jf(ca)  —  d||,  ||/?ca||)  lies  in  the  “corner”  of 
the  resulting  curve,  known  as  the  L- curve.  The  use  of  this  technique  for  determining  suitable 
choices  for  the  regularization  parameter  is  demonstrated  in  the  examples. 

In  each  of  the  three  examples,  the  data  was  sampled  on  a  regular  grid  {( xp ,  <,)}  C 
(0, 1)  x  (0,3].  Nine  equally  spaced  points  xp  =  pAx,Ax  =  .1,  were  taken  in  space  and  six 
equally  spaced  temporal  points  tq  =  qAt,  At  =  .5,  were  taken  for  a  total  of  n  =  54  data 
points.  To  the  data,  we  added  a  pseudo-random  noise  vector  e  from  a  Gaussian  distribution 
with  mean  0  and  standard  deviation  a  chosen  so  that  the  noise-to-signal  ratio  <7/||d||  =  .001; 
that  is,  noise  =  .1%  of  the  signal. 

The  summation  limits  were  taken  to  be  Af*  =  Ns  =  Mt  =  8  and  Nt  =  4  as  dictated  by 
the  choice  of  decay  parameters.  Hence  m,  =  17  basis  functions  were  used  in  the  expansion 
of  each  material  parameter.  The  startup  vectors  in  each  example  were  chosen  to  reflect  the 
positivity  and  general  endpoint  behavior  of  the  true  parameters.  All  problems  were  run  with 
sixteen  place  accuracy  on  a  Vax  8550. 
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Example  4-1- 

In  this  example,  the  parameter  El  is  considered  to  be  known  and  the  parameter  cpl 
is  to  be  numerically  recovered.  The  startup  vector  was  taken  to  be  the  mx  x  1  vector 
cq  =  [2,  .5,  .5,  •  •  ■ ,  .5,  .5,  —  2]T.  For  varying  values  of  the  regularization  parameter  a,  the 
L-curve  is  plotted  in  Figure  4.1.  Note  that  the  points  a  =  10-5  and  a  =  5  x  10-6  yield 
points  (||tf(ca)  —  d||,  || Rca || )  in  the  “corner”  of  the  curve.  For  a  =  5  x  10-4,  a  =  10~5  and 
a  =  10-7,  the  plots  of  the  true  and  approximate  damping  parameters  are  given  in  Figure  4.2. 
It  can  be  seen  that  the  “corner”  value,  a  =  10-s,  provides  a  good  choice  for  the  regularization 
parameter  whereas  a  =  10-7  is  not  large  enough  to  damp  out  the  contribution  due  to  the 
smaller  singular  values.  Finally,  the  choice  a  =  5  x  10”4  causes  too  much  smoothing  and 
information  about  the  parameter  is  lost.  The  results  from  this  example  demonstrate  the 
viability  of  the  method  for  problems  in  which  El  is  known  and  CqI  is  unknown. 


Figure  4.1:  The  Tikhonov  L-Curve  for  Example  4.1. 
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Figure  4.2:  True  and  Approximate  Damping  Parameters  for  Example  4.1 

- (a  =  5  x  10'4), - (a  =  lO'8),  •••(<*  =  10~7), - (True). 


Example  4-2 

Consider  now  the  case  where  the  stiffness  parameter  El  is  considered  to  be  unknown 
and  the  damping  parameter  coJ  is  assumed  to  be  known.  Here  the  mm  x  1  startup  vector 
Co  =  [1, 1,.5,  •  •  • ,  .5,2,0]t  was  used.  Since  the  L-curve  for  this  example  was  very  similar 
to  that  in  Example  4.1,  the  true  and  approximate  stiffness  parameters  corresponding  to 
a  =  5  x  10~4,  a  =  10"8  and  a  =  10-7  were  computed  with  results  given  in  Figure  4.3.  It  is 
again  noted  that  the  “corner”  value  of  a  =  10~8  provides  a  good  choice  for  the  regularization 
parameter  whereas  the  choice  a  =  5  x  10-4  causes  too  much  smoothing.  Finally,  the  error 
contributions  due  to  the  smaller  singular  values  become  more  apparent  with  a  =  10~7. 
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X-axis 


Figure  4.3:  True  and  Approximate  Stiffness  Parameters  for  Example  4.2 

- (a  =  5  x  lO-4), - (a  =  10~5),  •••(«  =  10“7), - (True). 

Example  4  $ 

In  this  example,  the  method  is  used  to  simultaneously  recover  both  the  stiffness  parameter 
El  and  the  damping  parameter  c©  J.  Following  from  the  previous  two  examples,  the  2m,  x  1 
initial  vector  had  the  form  Co  =  [c0,  co]T.  The  true  and  approximate  material  parameters 
corresponding  to  a  =  10~*  are  plotted  in  Figures  4.4  and  4.5.  From  the  figures  it  is  noted  that 
although  the  method  is  accurately  recovering  the  shape  of  the  functions,  there  is  more  error 
in  the  magnitude  than  in  the  previous  two  examples.  Hence  it  appears  that  a  larger  number 
of  state  and  parameter  basis  functions  are  needed  to  accurately  recover  both  parameters 
simultaneously,  and  current  efforts  are  directed  at  devising  linear  algebra  techniques  which 
would  facilitate  the  solution  of  the  larger  discrete  systems. 
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Figure  4.4:  True  and  Approximate  Stiffness  Parameters  for  Example  4.3 
- (a  =  10-*), - (True). 


Figure  4.5:  True  and  Approximate  Damping  Parameters  for  Example  4.3 
- (a  =  10~8), - (True). 
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